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Abstract 

An accurate simulation of the propagation of muons through matter is needed for the 
analysis of data produced by muon/neutrino underground experiments. A muon may sustain 
hundreds of interactions before it is detected by the experiment. Since a small systematic 
uncertainty repeated hundreds of times may lead to sizable errors, requirements on the 
precision of the muon propagation code are very stringent. A new tool for propagating 
muon and tau charged leptons through matter that is believed to meet these requirements is 
presented here. An overview of the program is given and some results of its application are 
discussed. 



1 Introduction 



In order to observe atmospheric and cosmic neutrinos with a large underground 
detector (e.g., AMANDA [1]), one needs to isolate the neutrino signal from the 
3-5 orders of magnitude larger signal from the background of atmospheric muons. 
Methods that do this have been designed and proven viable [2]. In order to prove 
that these methods work and to derive indirect results such as the spectral index 
of atmospheric muons, one needs to compare data to the results of the computer 
simulation. Such a simulation normally contains three parts: propagation of the 
measured flux of the cosmic particles from the top of the atmosphere down to the 
surface of the ground (ice, water); propagation of the atmospheric muons from the 
surface down to and through the detector; and generation of the secondary particles 
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(electrons, Cherenkov photons, etc.) in the vicinity of the detector and their inter- 
action with the detector components. The first part is normally called generator, 
since it generates muon flux at the ground surface; the second is propagator, and 
the third simulates the detector interaction with the passing muons. To generate at- 
mospheric muon and neutrino fluxes we used CORSIKA [3]. Results and methods 
of using CORSIKA as a generator in a neutrino detector (AMANDA-II) were dis- 
cussed in [4,5]. Several muon propagation Monte Carlo programs were used with 
different degrees of success as propagators. Some are not suited for applications 
which require the code to propagate muons in a large energy range (e.g., mudedx, 
a.k.a. LOH [6]), and the others seem to work in only some of the interesting energy 
range (E > 1 TeV, propmu, a.k.a. LIP [7]) [8]. Most of the programs use cross 
section formulae whose precision has been improved since their time of writing. 
For some applications, one would also like to use the code for the propagation of 
muons that contain 100 — 1000 interactions along their track, so the precision of 
each step should be sufficiently high and the computational errors should accumu- 
late as slowly as possible. Significant discrepancies between the muon propagation 
codes tested in this work were observed, and are believed to be mostly due to algo- 
rithm errors (see Appendix B). This motivated writing of a new computer program 
(Muon Monte Carlo: MMC [9]), which minimizes calculational errors, leaving only 
those uncertainties that come from the imperfect knowledge of the cross sections. 



2 Description of the code 

The primary design goals of MMC were computational precision and code clarity. 
The program is written in Java, its object-oriented structure being used to improve 
code readability. MMC consists of pieces of code (classes), each contained in a sep- 
arate file. These pieces fulfill their separate tasks and are combined in a structured 
way (Figure 1). 

The code evaluates many cross-section integrals, as well as several tracking inte- 
grals. All integral evaluations are done by the Romberg method of the 5th order (by 
default) [10] with a variable substitution (mostly log-exp). If an upper limit of an 
integral is an unknown (that depends on a random number), an approximation to 
that limit is found during normalization integral evaluation, and then refined by the 
Newton-Raphson method combined with bisection [10]. 

Originally, the program was designed to be used in the Massively Parallel Network 
Computing (SYMPHONY) [11] framework, and therefore computational speed 
was considered only a secondary issue. However, parametrization and interpolation 
routines were implemented for all integrals. These are both polynomial and rational 
function interpolation routines spanned over a varying number of points (5 by de- 
fault) [10]. Inverse interpolation is implemented for root finding (i.e., when x(f) is 
interpolated to solve f(x) = y). Two-dimensional interpolations are implemented 
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Fig. 1. MMC structure 



as two consecutive one-dimensional ones. It is possible to turn parameterizations 
on or off for each integral separately at program initialization. The default energy 
range in which parametrized formulae will work was chosen to be from 105.7 MeV 
(the muon rest mass; 1777 MeV for taus) to Eu g = 10 14 MeV, and the program was 
tested to work with much higher settings of E biy . With full optimization (parameter- 
izations) this code is at least as fast or even faster than the other muon propagation 
codes discussed in Appendix B. 

Generally, as a muon travels through matter, it loses energy due to ionization losses, 
bremsstrahlung, photo-nuclear interaction, and pair production. The cross section 
formulae are summarized in Section 9. These formulae are claimed to be valid to 
within about 1% in the energy range up to > 10 TeV Theoretical uncertainties in 
the photonuclear cross section above 100 TeV are higher. All of the energy losses 
have continuous and stochastic components, the division between which is artifi- 
cial and is chosen in the program by selecting an energy cut (e cut , also E cut ) or a 
relative energy loss cut (v cut ). In the following, v cut and e cut are considered to be 
interchangable and related by e cut = v cut E (even though only one of them is a con- 
stant). Ideally, all losses should be treated stochastically. However, that would bring 
the number of separate energy loss events to a very large value, since the probabil- 
ity of such events to occur diverges as 1/ E lost for the bremsstrahlung losses, as 
the lost energy approaches zero, and even faster than that for the other losses. In 
fact, the reason this number, while being very large, is not infinite, is the existence 
of kinematic cutoffs (larger than some e ) for all diverging cross sections. A good 
choice of v cut for the propagation of atmospheric muons should lie in the range 
0.05 — 0.1 (Section 3, also [12]). For monoenergetic beams of muons, v cut may 
have to be chosen to be high as 10~ 3 — 10~ 4 . 
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2.1 Tracking formulae 



Let the continuous part of the energy losses (a sum of all energy losses, integrated 
from zero to e cut ) be described by a function f(E): 

-§-'<■>■ 



X; dx X, 

Fig. 2. Derivation of tracking formulae 

The stochastic part of the losses is described by the function a(E), which is a prob- 
ability for any energy loss event (with lost energy > e cut ) to occur along a path of 1 
cm. Consider the particle path from one interaction to the next consisting of small 
intervals (Figure 2). On each of these small intervals the probability of interaction 
is dP(E(xi)) = a(E(xi))dx. We now derive an expression for the final energy 
after this step as a function of the random number £. The probability to completely 
avoid stochastic processes on an interval (Xi,Xf) and then suffer a catastrophic loss 
on dx at xj is 

(1 - dP(E( Xi ))) ■ ... • (1 - dP{E{x f - X ))) ■ dP(E(x f )) 
« exp(-dP(E(xi))) ■ ... • exp(-dP(E(x f -i))) ■ dP(E(x f )) 
exp (- Q dP(E(x))^ ■ dP(E(x f )) 

= d f f-exp(- Q • dE)\ = d(-0, e G (0; 1] 

To find the final energy after each step the above equation is solved for Ef. 



f E f a(E) 

/ - dE = - log(0 (energy integral). 



This equation has a solution if 



(: I o i 



f(E) 



Here ti ow is a low energy cutoff, below which the muon is considered to be lost. 
Note that f(E) is always positive due to ionization losses (unless e cut < 
The value of a(E) is also always positive because it includes the positive decay 
probability. If £ < £ , the particle is stopped and its energy is set to ei ow . The 
corresponding displacement for all £ can be found from 

r E f dE 

x f = Xi- / -jrr— (tracking integral), 

•ie, f{E) 
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and time elapsed can be found from 

r x f dx r E f dE 

tf = U+ / -j-^ =U- (time integral). 

Jxi v(x) Je, f(E)v(E) 

Evaluation of time integral based on the approximation v = c,tf = ti + (xj — Xi)/c, 
is also possible. 



2.2 Continuous randomization 



It was found that for higher v cut muon spectra are not continuous (Figure 3). In 
fact, there is a large peak (at E pea k) that collects all particles that did not suffer 
stochastic losses followed by the main spectrum distribution separated from the 
peak by at least the value of v cut E pea k (the smallest stochastic loss). The appearance 
of the peak and its prominence are governed by v cut , co-relation of initial energy 
and propagation distance, and the binning of the final energy spectrum histogram. 
In order to be able to approximate the real spectra with even a large v cut and to 
study the systematic effect at a large v cut , a continuous randomization feature was 
introduced. 
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Fig. 3. Distribution of the final energy 
of the muons that crossed 300 m of 
Frejus Rock with initial energy 100 TeV: 
Vcut = 0.05 (solid), v cut = 10" 4 
(dashed-dotted), v cu t = 0.05 and cont 
option (dotted) 



Fig. 4. A close-up on the Figure 3: 
Vcut = 0.05 (solid), Vcut = 0.01 
(dashed), v cut = 10~ 3 (dotted), 
Vcut = 10~ 4 (dashed-dotted) 



For a fixed v cut or e cu t a particle is propagated until the algorithm discussed above 
finds an interaction point, i.e., a point where the particle loses more than the cutoff 
energy. The average value of the energy decrease due to continuous energy losses 
is evaluated according to the energy integral formula of the previous section. There 
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Fig. 5. Same as in Figure 4, but with 
cont option enabled 



Fig. 6. Comparison of parametrized 
(dashed-dotted) with exact (non- 
parametrized, dotted) versions for 
v C ut = 0.01. Also shown is the relative 
difference of the curves, i.e., ratio of 
the difference over exact histogram bin 
values. 



will be some fluctuations in this energy loss, which are not described by this for- 
mula. Let us assume there is a cutoff for all processes at some small e <C e cut . 
Then the probability p(e; E) for a process with e < ei ost < e cut on the distance dx 
is finite. Now choose dx so small that 

re c nt 

Po — p( e ', E) fie • dx <C 1. 

Jeo 

Then the probability to not have any losses is 1 — p , and the probability to have 
two or more separate losses is negligible. The standard deviation of the energy loss 
on dx from the average value 

re cu t 

< e >= / e ■ p(e; E) de ■ dx 

J e 

is then < (Ae) 2 >=< e 2 > — < e > 2 , where 

r e cut 

<e >= / e • p(e; E) de ■ dx. 

Je 



If the value of v cut or e cut used for the calculation is sufficiently small, the distance 
xj — Xi determined by the energy and tracking integrals is so small that the average 
energy loss E { — Ef is also small (as compared to the initial energy Ei). One may 
therefore assume p(e; E) ~ p(e; Ej), i.e., the energy loss distributions on the small 
intervals dx n that sum up to the Xf — Xi, is the same for all intervals. Since the 
total energy loss Ei — Ef = e„, the central limit theorem can be applied, and the 
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final energy loss distribution will be Gaussian with the average AE = Ei — Ej and 
width 

< (A(A£)) 2 >= £ (< el > - < e n > 2 ^ 



= E 



e 2 n ■ p{e n ; E { ) de n ^j dx n - Qf e„ • p(e n ; c?e„^) dx 2 



~ / dx • i e • p(e; E{x)) de) — dx • I / e • p(e; E{x)) de ) dx 

JXi \Jeo / JXi \Jeo / 

Here was replaced with the average expectation value of energy at E'(x). As 
dx — > 0, the second term disappears. The lower limit of the integral over e can be 
replaced with zero, since none of the cross sections diverge faster than or as fast as 
1/e 3 . Then, 

r x f dE / f e cut \ 
< (A(AE)) >~ J _^ • ^ e 2 -p(e;E) de) (co/rf integral). 

This formula is applicable for small v cut , as seen from the derivation. Energy spectra 
calculated with continuous randomization converge faster than those without as v cut 
is lowered (see Figures 4 and 5). 



3 Computational and algorithm errors 



All cross-section integrals are evaluated to the relative precision of 10~ 6 ; the track- 
ing integrals are functions of these, so their precision was set to a larger value of 
1CT 5 . To check the precision of interpolation routines, results of running with pa- 
rameterizations enabled were compared to those with parameterizations disabled. 
Figure 7 shows relative energy losses for ice due to different mechanisms. De- 
cay energy loss is shown here for comparison and is evaluated by multiplying the 
probability of decay by the energy of the particle. In the region below 1 GeV, 
bremsstrahlung energy loss has a double cutoff structure. This is due to a differ- 
ence in the kinematic restrictions for muon interaction with oxygen and hydrogen 
atoms. A cutoff (for any process) is a complicated structure to parametrize and 
with only a few parametrization grid points in the cutoff region, interpolation er- 
rors (e pa — e np ) I e pa may become quite high, reaching 100% right below the cutoff, 
where the interpolation routines give non-zero values, whereas the exact values are 
zero. But since the energy losses due to either bremsstrahlung, photonuclear pro- 
cess, or pair production are very small near the cutoff in comparison to the sum of 
all losses (mostly ionization energy loss), this large relative error results in a much 
smaller increase of the relative error of the total energy losses (Figure 8). Because 
of that, parametrization errors never exceed 1CT 4 — 1CT 3 , for the most part being 
even much smaller (1CT 6 — 1CT 5 ), as one can estimate from the plot. These errors 
are much smaller than the uncertainties in the formulae for the cross sections. Now 
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the question arises whether this precision is sufficient to propagate muons with hun- 
dreds of interactions along their way. Figure 6 is one of the examples that demon- 
strate that it is sufficient: the final energy distribution did not change after enabling 
parametrizations. Moreover, different orders of the interpolation algorithm (g, cor- 
responding to the number of the grid points over which interpolation is done) were 
tested (Figure 9) and results of propagation with different g compared with each 
other (Figure 10). The default value of g was chosen to be 5, but can be changed to 
other acceptable values 3 <g< 6 at the run time. 
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Fig. 7. Ionization (upper solid curve), 
bremsstrahlung (dashed), photonuclear 
(dotted), epair production (dashed-dot- 
ted) and decay (lower solid curve) losses 
in ice 
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Fig. 9. Interpolation precision for differ- 
ent orders of the interpolation algorithm 



Fig. 10. Comparison of the result of the 
propagation for different orders of the in- 
terpolation algorithm 



MMC employs a low energy cutoff e [ow below which the muon is considered to be 
lost. By default it is equal to the mass of the muon, but can be changed to any higher 
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value. This cutoff enters the calculation in several places, most notably in the initial 
evaluation of the energy integral. To determine the random number £ below which 
the particle is considered stopped, the energy integral is first evaluated from E{ to 
e\ ow . It is also used in the parametrization of the energy and tracking integrals, since 
they are evaluated from this value to Ei and Ef, and then the interpolated value for 
Ef is subtracted from that for E{. Figure 11 demonstrates the independence of 



MMC from the value of ei ow . For the curve with t\ c 



m, 



integrals are evaluated 



in the range 105.7 MeV — 100 TeV, i.e., over six orders of magnitude, and they are 
as precise as those calculated for the curve with ei ow =\0 TeV, with integrals being 
evaluated over only one order of magnitude. 
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Fig. 11. Comparison of ei ow = (dot- 
ted-dashed) with e/ otu =10 TeV (dotted). 
Also shown is the relative difference of 
the curves. 
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Fig. 12. Ionization (upper solid curve), 
bremsstrahlung (dashed), photonuclear 
(dotted), epair production (dashed-dot- 
ted) spectra for £^=10 TeV in the Frejus 
rock 



Figure 12 demonstrates the spectra of secondaries (delta electrons, bremsstrahlung 
photons, excited nuclei, and electron pairs) produced by the muon, whose energy 
is kept constant at 10 TeV. The thin lines superimposed on the histograms are the 
probability functions (cross sections) used in the calculation. They have been cor- 
rected to fit the logarithmically binned histograms (multiplied by the size of the 
bin which is proportional to the abscissa, i.e., the energy). While the agreement is 
trivial from the Monte Carlo point of view, it demonstrates that the computational 
algorithm is correct. 

Figure 13 shows the relative deviation of the average final energyof the 4-10 6 1 TeV 
and 100 TeV muons propagated through 100 m of Frejus RocliU with the abscissa 
setting for v cut , from the final energy obtained with v cut = 1. lust like in [12] the 
distance was chosen small enough so that only a negligible number of muons stop, 



1 A medium with properties similar to that of standard rock (see second table in Appendix 
A) used for data analysis in the Frejus experiment [13]. 
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while large enough so that the muon suffers a large number of stochastic losses 
(> 10 for v cut < 10~ 3 ). All points should agree with the result for v cut = 1, since it 
should be equal to the integral of all energy losses, and averaging over the energy 
losses for v cut < 1 is evaluating such an integral with the Monte Carlo method. 
There is a visible systematic shift < (1 — 2) ■ 10 -4 (similar for other muon energies), 
which can be considered as another measure of the algorithm accuracy [12]. 




Fig. 13. Algorithm errors (average energy loss) 
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Fig. 14. 10 6 muons with energy 9 TeV propagated through 10 km of water: regular 
(dashed) vs. cont (dotted) 



In the case when almost all muons stop before passing the requested distance (see 
Figure 14), even small algorithm errors may substantially affect survival proba- 
bilities. Table 1 summarizes the survival probabilities for a monochromatic muon 
beam of 10 6 muons with three initial energies (1 TeV, 9 TeV, and 10 6 TeV) going 
through three distances (3 km, 10 km, and 40 km) in water. One should note that 
these numbers are very sensitive to the cross sections used in the calculation; e.g., 
for 10 9 GeV muons propagating through 40 km the rate increases 12% when the 
BB1981 photonuclear cross section is replaced with the ZEUS parametrization (see 
Figure 37). However, the same set of formulae was used throughout this calcula- 
tion. The errors of the values in the table are statistical and are < ±0.001. The 
survival probabilities converge on the final value for v cut < 0.01 in the first two 
columns. Using the cont option helped the convergence in the first column. How- 
ever, the cont values departed from regular values more in the third column. The 
relative deviation (5.4%) can be used as an estimate of the continuous randomiza- 
tion algorithm precision (not calculational errors) in this case. One should note, 
however, that with the number of interactions > 10 3 the continuous randomization 
approximation formula was applied > 10 3 times. It explains why the value of cont 
version for v cut = 0.01 is closer to the converged value of the regular version than 
for v cut = 10~ 3 . 
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Table 1 

Survival probabilities 



Vcut 


cont 


1 TeV 3 km 


9 TeV 10 km 


10 6 TeV 40 km 


0.2 


no 








0.081 


0.2 


yes 


0.009 


0.052 


0.113 


0.05 


no 





0.028 


0.076 


0.05 


yes 


0.041 


0.034 


0.073 


0.01 


no 


0.027 


0.030 


0.075 


0.01 


yes 


0.031 


0.030 


0.072 


10~ 3 


no 


0.031 


0.031 


0.074 


10~ 3 


yes 


0.031 


0.030 


0.070 



4 Electron, tau, and monopole propagation 
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Fig. 16. Tau energy losses in Frejus Rock 

Electrons and taus can also be propagated with MMC. Bremsstrahlung is the domi- 
nant cross section in case of electron propagation, and the complete screening case 
cross section should be selected (Section 9.2.4). Electron energy losses in Ice are 
shown in Figure 15 (also showing the LPM suppression of cross sections). 

For tau propagation Bezrukov-Bugaev parameterization with the hard component 
(Section 9.3.1) or the ALLM parametrization (Section 9.3.2) should be selected 
for photonuclear cross section. Tau propagation is quite different from muon prop- 
agation because the tau lifetime is 7 orders of magnitude shorter than the muon 
lifetime. While muon decay can be neglected in most cases of muon propaga- 
tion, it is the main process to be accounted for in the tau propagation. Figures 
16 and 17 compare tau energy losses with losses caused by tau decay (given by 









ioniz 


























brcm 
phou 


) 


























total 
































































































































































































• 



























: ^ ^ Y _ _ .i _ _ 

-3 -2 -1 2345678910 

10 10 TO 1 10 10 10 10 10 10 10 10 1010 



energy [ GeV ] 

Fig. 15. Electron energy losses in Ice 
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Fig. 18. Average range of taus propa- 
gated through Frejus Rock 

E T J (pv T r) = m T j (pv T r ); this is the energy per mwe deposited by decaying taus 
in a beam propagating though a medium with density p). Figure 18 compares the 
average range of taus propagated through Frejus Rock with v cut = 1 (completely 
continuously) and v cut = 1CT 3 (detailed stochastic treatment). Both treatments pro- 
duce almost identical results. Therefore, tau propagation can be treated continu- 
ously for all energies unless one needs to obtain spectra of the secondaries created 
along the tau track. 

Monopoles can also be propagated with MMC. All cross sections except brems- 
strahlung (which scales as z 4 ) are scaled up with a factor z 2 , where z = I /(2a) is 
the monopole charge, according to [14,15]. 



5 Comparison with other propagation codes 

Several propagation codes have been compared with MMC. Where possible MMC 
settings were changed to match those of the other codes. Figure 20 compares energy 
losses calculated with MMC and MUM [12], and Figure 19 compares the results of 
muon propagation through 800 m of ice with MMC and MUM (v cut = 10~ 3 , ZEUS 
parametrization of the photonuclear cross section, Andreev Berzrukov Bugaev pa- 
rameterization of bremsstrahlung). 

Survival probabilities of Table 1 were compared with results from [12] in Table 2. 
Survival probabilities are strongly correlated with the distribution of the highest- 
energy muons in an originally monoenergetic beam. This, in turn, is very sensitive 
to the algorithm errors and the cross-section implementation used for the calcula- 
tion. 

A detailed comparison between spectra of secondaries produced with MMC, MUM, 
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Fig. 17. Sum of tau energy losses in 
Frejus Rock 
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Fig. 19. Comparison of energy losses in 
water. Red MMC curve shows photonu- 
clear losses with v integration limits as 
prescribed in [12] 



Fig. 20. Comparison of muon propaga- 
tion through 800 m of Ice with MMC and 
MUM 



Table 2 

Survival probabilities of MMC compared to other codes 



Vcut 


propagation code 


1 TeV 3 km 


9 TeV 10 km 


10 6 TeV 40 km 


io- 3 


MMC (BB81) 


0.031 


0.031 


0.074 


io- 3 


MMC (ZEUS) 


0.031 


0.030 


0.083 


10~ 3 


MUM [12] 


0.029 


0.030 


0.078 


10~ 3 


MUSIC [16] 


0.033 


0.031 


0.084 


10~ 3 


PROPMU [17] 


0.19 


0.048 


0.044 



LOH [6], and LIP [7,17] is given in the Appendix B. A definite improvement of 
MMC over the other codes can be seen in the precision of description of spectra of 
secondaries and the range of energies over which the program works. 



6 Energy losses in ice and rock, some general results 



The code was incorporated into the Monte Carlo chains of three detectors: Frejus 
[13,18], AMANDA [8,4], and IceCube [19]. In this section some general results 
are presented. 
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6.1 Average muon energy loss 



The plot of energy losses was fitted to the function dE/dx = a + bE (Figure 21). 
The first two formulae for the photonuclear cross section (Section 9.3.1) can be 







dH/dx=i 
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Fig. 21. Fit to the energy losses in ice Fig. 22. x 2 plot for energy losses in ice 

fitted the best, all others lead to energy losses deviating more at higher energies 
from this simple linear formula; therefore the numbers given were evaluated using 
the first photonuclear cross section formula. In order to choose low and high energy 
limits correctly (to cover the maximum possible range of energies that could be 
comfortably fitted with a line), a x 2 plot was generated and analyzed (Figure 22). 
The green curve corresponds to the x 2 °f the fit with a fixed upper bound and a 
varying lower bound on the fitted energy range. Correspondingly, the blue curve 
describes the x 2 °f the fit with a fixed lower bound and a varying upper bound. 
The x 2 a t low energies goes down sharply, then plateaus at around 10 GeV. This 
corresponds to the point where the linear approximation starts to work. For the high 
energy boundaries, x 2 rises monotonically. This means that a linear approximation, 
though valid, has to describe a growing energy range. An interval of energies from 
20 GeV to 10 11 GeV is chosen for the fit. Table 3 summarizes the found fits to 
a and b; the errors in the evaluation of a and b are in the last digit of the given 
number. However, if the lower energy boundary of the fitted region is raised and/or 
the upper energy boundary is lowered, each by an order of magnitude, a and b 
change by about 1%. 

To investigate the effect of stochastic processes, muons with energies 105.7 MeV 
— 10 11 GeV were propagated to the point of their disappearance. The value of 
v C ut — 5 • 10~ 3 was used in this calculation; using the the continuous randomization 
option did not change the final numbers. The average final distance (range) for each 
energy was fitted to the solution of the energy loss equation dE/dx = a + bE: 

x f = log(l + Ei ■ b/a)/b 
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Table 3 



Fits to a and b for continuous losses (average energy losses) 



medium 


a GeV 
mwe 


b 10-3 
' mwe 


av. dev. 


max. dev. 


a GeV 
mwe 


b 10-3 
' mwe 


a GeV 
mwe 


b 10-3 
mwe 




20 - 10 11 GeV 


20 - 10 7 GeV 


ALLM97 


air 


0.281 


0.347 


3.6% 


6.5% 


0.284 


0.335 


0.282 


0.344 


ice 


0.259 


0.363 


3.7% 


6.6% 


0.262 


0.350 


0.260 


0.360 


fr. rock 


0.231 


0.436 


3.0% 


5.1% 


0.233 


0.423 


0.231 


0.431 


st. rock 


0.223 


0.463 


2.9% 


5.1% 


0.225 


0.451 


0.224 


0.459 



(Figure 23). The same analysis of the x 2 plot as above was done in this case (Figure 
24). A region of initial energies from 20 GeV to 10 11 GeV was chosen for the fit. 
Table 4 summarizes the results of these fits. 

Table 4 

Fits to a and b for stochastic losses (average range estimation) 



medium 


a GeV 
' mwe 


b 10-3 
mwe 


av. dev. 


ice 


0.268 


0.470 


3.0% 


frejus rock 


0.218 


0.520 


2.8% 
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Fig. 23. Fit to the average range in Frejus 
rock 
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Fig. 24. x 2 pl ot f° r average range in 
Frejus rock 



As the energy of the muon increases, it suffers more stochastic losses before it is 
los{3 and the range distribution becomes more Gaussian-like (Figure 31). It is also 
shown in the figure (vertical lines) that the inclusion of stochastic processes makes 
the muons on average travel a shorter distance. 



As considered by the algorithm, here: stopped. 
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6.2 Muon range 



In certain cases it is necessary to find the maximum range x of (the majority of) 
muons of certain energy E, or find what is the minimum energy E cut muons must 
have in order to cross distance x. 




distance z [ km ] fraction of muons left after passing distance z 



r . » c „ «• . —i ,. • Fie. 26. Distance in ice vs. fraction of 

Fig. 25. Muon range distributions m ice fc 

survived muons 

To determine such function E cut (x), MMC was run for ice as propagation medium, 
with muon energies from 105 MeV to 10 20 eV. For each energy 10 5 muons were 
propagated to the point of their disappearance and the distance traveled was his- 
togrammed (Figure 25). This is similar to the analysis done in Section 6.1. How- 
ever, instead of the average distance traveled, the distance at which only a frac- 
tion of muons survives was determined for each muon energy (Figure 26). Two 
fixed fractions were used: 99% and 99.9%. MMC was run with 2 different settings: 
v cut = 1CT 2 with the cont {continuous randomization feature described in Section 
2.2) option and v cut = without cont. In Figure 27 the ratio of distances de- 
termined with both settings is displayed for 99% of surviving muons (red line) and 
for 99.9% (green line). Both lines are very close to 1.0 in most of the energy range 
except the very low energy part (below 2 GeV) where the muon does not suffer 
enough interactions with the v cut = 1CT 2 setting before stopping (which means 
v cu t has to be lowered for a reliable estimation of the shape of the travelled distance 
histogram). The ratio of 99% distance to 99.9% distance is also plotted (dark and 
light blue lines). This ratio is within 10% of 1, i.e., 0.1% of muons travel less than 
10% farther than 1% of muons. 

The value v cut = 10~ 3 with no cont setting, used to determine the maximum range 
of the 99.9% of the muons, was chosen for the estimate of the function E cut (x). 
The function 

x f = \og(l+Ei-b/a)/b, 
which is a solution to the equation represented by the usual approximation to the 
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muon energy [ GeV ] fit range boundary [ GeV ] 



Fig. 27. Comparison between different 

surviving fraction and MMC configura- Fi §" 28 " * of the fit as a function of fit 
tion settings boundaries 

energy losses: dE/dx = a + bE, was fitted to E cut (x). Figure 28 shows the x 2 
of the fit as function of the lower (green) and upper (blue) boundaries of the fitted 
energy range. Using the same argument as in Section 6.1 the lower limit is chosen 
at just below 1 GeV while the upper limit was left at 10 11 GeV. As seen from the 
plot, raising the lower boundary to as high as 400 GeV would not lower the x 2 of 
the fit (and the root mean square of the deviation from it), so the lower boundary 
was left at 1 GeV for generality of the result. The fit is displayed in Figure 29 and 
the deviation of the actual Xf from the fit is shown on Figure 30. The maximum 
deviation is less than 20%, which can be accounted for by lowering a and b by 
20%. Therefore, the final values quoted here for the function 

E cut (x) = (e bx - l)a/b 

GeV , , , 1 

are a = 0.212/1.2 and b = 0.251 ■ lCT 3 / 1 - 2 



mwe mwe 
The distances obtained with these values for four different muon energies are shown 
by red solid lines in Figure 25. The distances obtained with values of a and b not 
containing the 20% correction are shown with green dashed lines. 



7 Phenomenological lepton generation and neutrino propagation 



MMC allows one to generate fluxes of atmospheric leptons according to parame- 
terizations given in [20]. Earth surface (important for detectors at depth) and atmo- 
spheric curvature are accounted for, and so are muon energy losses and probability 
of decay. Although the reference [20] provides flux parameterization, which is ac- 
curate in the region of energies from 600 GeV to 60 TeV, it is possible to introduce 
a correction to spectral index and normalization of each leptonic component and 
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dE/dx=a+bE 
a= 0.212 [GeV/mwe], 
b= 0.251 • 10" 3 [ 1/mwe] 
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Fig. 29. Fittothe£ CMt (:r) 



Fig. 30. Deviation of the E cut (x) from 
the fit 
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Fig. 32. 3 regions of propagation defined 
for AMANDA-II simulation 



Fig. 31. Range distributions in Frejus 
rock: solid line designates the value of 
the range evaluated with the second table 
(continuous and stochastic losses) and 
the broken line shows the range eval- 
uated with the first table (continuous 
losses only). 



extrapolate the results to the desired energy range. One can also add an ad-hoc 
prompt component, specify _E ,_7 -like fluxes of neutrinos of all flavors, or inject 
leptons with specified location and momenta into the simulation. 



Neutrino cross sections are evaluated according to [21,22,23] with CTEQ6 par- 
ton distribution functions [24] (Figure 33). Neutrino and anti-neutrino neutral and 
charged current interaction, as well as Glashow resonance z/ e e~ cross sections are 
taken into account. Power-law extrapolation of the CTEQ PDFs to small x is im- 
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plemented to extend the cross section applicability range to high energies. Earth 
density is calculated according to [25,22], with a possibility of adding layers of dif- 
ferent media. All secondary leptons are propagated, therefore it is possible to simu- 
late particle oscillations, e.g., r <-> v T . Additionally, atmospheric neutrino v T <-> 
oscillations are simulated (Figure 34). 
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Fig. 33. Simulated neutrino cross sec- 
tions: higher blue curves are CC, lower 
red curves are NC; solid are v, dashed 
are v; green dotted is v e e~ — > W~ 



Fig. 34. Neutrino flavor oscillations: 
muon events. Since latitude-dependent 
geomagnetic cutoff is not calculated, a 
fixed 10 GeV cutoff is applied (cf. [26]). 



Alternatively, MMC allows integration with other neutrino generators/propagators, 
such as NUSIM [27], ANIS [26], or Juliet [28], or lepton generators, e.g., COR- 
SIKA [3]. 



8 MMC implementation for AMANDA-II 



Most light observed by AMANDA-II is produced by muons passing through a 
cylinder with radius 400 and length 800 meters around the detector [4]. Inside 
this cylinder, the Cherenkov radiation from the muon and all secondary showers 
along its track with energies below 500 MeV (a somewhat loose convention) are 
estimated together. In addition to light produced by such a "dressed" muon, all sec- 
ondary showers with energies above 500 MeV produced in the cylinder create their 
own Cherenkov radiation, which is considered separately for each secondary. So 
in the active region of the detector muons are propagated with E cut = 500 MeV, 
creating secondaries on the way. This is shown as region 2 in the Figure 32. 

In region 1, which is where the muon is propagated from the Earth's surface (or 
from under the detector) to the point of intersection of its track with the detector 
cylinder, muons should be propagated as fast as possible with the best accuracy. 
For downgoing muons, values of v cut = 0.05 with the continuous randomization 
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option enabled were found to work best. These values should also work for muons 
propagated from points which are sufficiently far from the detector. For muons 
created in the vicinity of the detector, values of v cut = 0.01 with cont or even 
v C ut = 0.001 without cont should be used. 

In region 3, which is where the muon exits the detector cylinder, it is propagated in 
one step (v cut = 1.0, no cont) to the point of its disappearance, thus only resulting 
in an estimate of its average range. 

It is possible to define multiple concentric media to describe both ice and rock be- 
low the ice, which is important for the study of the muons, which might be created 
in either medium in or around the detector and then propagated toward it. Defini- 
tion of spherical, cylindrical, and cuboid detector and media geometries is possible. 
This can be easily extended to describe other shapes. 

Although the ALLM97 with nuclear structure function as described in Section 9.3.3 
parametrization of the photonuclear cross section was chosen to be the default for 
the simulation of AMANDA-II, other cross sections were also tested. No signif- 
icant changes in the overall simulated data rate or the number of channels (N ch ) 
distribution (important for the background muon analysis of [4,29]) were found be- 
tween the parameterizations described in Section 9.3. This is to be expected since 
for the background muons (most of which have energies of 0.5-10 TeV on the sur- 
face) all photonuclear cross section parameterizations are very close to each other 
(see Figure 37). Also the effects of the Moliere scattering and LPM-related effects 
(Section 9.5) can be completely ignored (although they have been left on for the 
default settings of the simulation). 



9 Formulae 



This section summarizes cross-section formulae used in MMC. In the formulae 
below, E is the energy of the incident muon, while v = vE is the energy of the 
secondary particle: knock-on electron for ionization, photon for bremsstrahlung, 
virtual photon for photonuclear process, and electron pair for the pair production. 
As usual, (3 = v jc and 7 = (1 — (3 2 ) ; also jj, is muon mass (or tau mass, except 
in the expression for q c of Section 9.2.3, where fi is just a mass-dimension scale 
factor equal to the muon mass [30]), m = m e is electron mass, and M is proton 
mass. Values of constants used below are summarized in Appendix A. 
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9.1 Ionization 



A standard Bethe-Bloch equation [31] was modified for muon and tau charged 
leptons (massive particles with spin 1/2 different from electron) following the pro- 
cedure outlined in [32]. The result is given below (and is consistent with [33]): 



dE 

dx 



Kz' 



Af3 2 



i(zy 



i 



upper 



where v„ 



2 V2£(l + l/7), 
2m e ( 7 2 - 1) 



l + 2 7 ^ 

The density correction 5 is computed as follows: 
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1 + 
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5 = 5 1 2{x ~ Xo) , ifX<X 

5 = 2(\nlO)X + C + a(X 1 - X) m , if X < X < X 1 

<5 = 2(lnlO)X + C, ifX>X 1? where X = log 10 (/3 7 ) 
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dudx 



l -Kz^^ 

2 A f3 2 v 2 



1 
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+ 



V 



2 V£(l + l/7), 



This formula, integrated from v min = • {^jff) to v U pper, gives the expression 
for energy loss above, less the density correction and (3 2 terms (plus two more terms 
which vanish if v min < u upper ). 



9.2 Bremsstrahlung 



According to [34], the bremsstrahlung cross section may be represented by the sum 
of an elastic component (a e i, discussed in [35,36]) and two inelastic components 

(A<), 

a = a d + Ao? + A<". 
9. 2. 1 Elastic Bremsstrahlung ( Kelner Kokoulin Petrukhin parameterization ): 
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a el (E,v) = -\ 2Z-r e 



4 4 2 
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3 3 
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2 a n 



where 5 ~ — 
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is the minimum momentum transfer. The formfactors (atomic A„ and nuclear A^) 
are 



A? (*) = In 
A?(*)=ln 



1 + 



Sy/eBZ-W/m, 
D„ 



l + 5(D n ^-2)/pi_ 
Integration limits for this cross section are 



LL = 1.54 A 1 



0.27 



< v < v mnx = 1 - hll^z^ 

4 # 



9.2.2 Petrukhin Shestakov form factor parameterization: 

Somewhat older parameterization of the form factors in the Bethe-Heitler formula 
[35] is given in [37]: 



a e i(E,v) = - 2Z—r e 

v \ fi 



4 4 
3 ~ 3 



v + v 2 ) (f)(5), 



(f)(5) = In 
(f)(5) = In 



189^ £-1/3 



1 + 189^^-1/3 
m 

2 189^ £-2/3 

3 m 



, £ < io 

, Z > 10. 



1 + I^^-l/S 
m 

9.2.3 Andre ev Berzrukov Bugaev parameterization: 

Another parameterization of the bremsstrahlung cross section, both elastic and in- 
elastic //-diagram contributions (not the e-diagram, which is included with the ion- 
ization cross section) is implemented according to [38,39,12]. 



a(E,v) = a \2r e Z 



m e \ 1 
H ) v 



[2-2v + w 2 )^i(g™„, Z) - -(1 - v)^ 2 (q mm , Z) 
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^1,2 (9min, Z) — ^i 2 (q m in, Z) — Ai >2 (q min , Z) 



*?(g min ,Z) = - 1 + ln 



H 2 aj 
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— #1 arctan h — 
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« 2 /\ 1 3, X( 
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9.2.4 Complete screening case: 

This parameterization is given in [31] (based on [40,41]) and is most suitable for 
electrons: 



a el {E,v) = ^[2Zy- e 



v + v' 

3 3 



Z 2 (L rad -f(Z)) + ZL' r 



rad 



+ 



1(1 - v)(Z 2 + Z)} , 



f(Z) = a 2 



1 + a 2 



+ 0.20206 - 0.0369a 2 + 0.0083a 4 - 0.002a 6 



, with a = ctZ. 



All bremsstrahlung parameterizations are compared in Figures 35 and 36. Param- 
eterization of Section 9.2.3 (abb) agrees best with the complete screening case of 
electrons and with the other two cross sections for muons, thereby providing the 
most comprehensive description of bremsstrahlung cross section. 
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Fig. 35. Bremsstrahlung cross section 
parameterizations for muons 
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Fig. 36. Bremsstrahlung cross section 
parameterizations for electrons 



9.2.5 Inelastic Bremsstrahlung: 



The effect of nucleus excitation can be evaluated as 



Bremsstrahlung on the atomic electrons can be described by the diagrams below; 
e-diagram is included with ionization losses (because of its sharp 1/v 2 energy loss 
spectrum), as described in [42]: 



d 2 N 



f d 2 N \ 



dvdx \ dvdx J 



Io 



a 
2^ 



(a(2b + c) - b 2 ) 



a = log(l + 2v/m e 



log((l - u/u max )/(l - u/E)), 



log((2 7 (l - v/E)m e )/{m„v/E)). 



The maximum energy lost by a muon is the same as in the pure ionization (knock- 
on) energy losses. The minimum energy is taken as is min = I(Z). In the above 
formula v is the energy lost by the muon, i.e., the sum of energies transferred to 
both electron and photon. On the output all of this energy is assigned to the electron. 
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Mi 



fi diagrams 



e diagrams 

The contribution of the yU-diagram (included with bremsstrahlung) is discussed in 
[34]: 

Aal n (E,v)^^-(2Zy e ) (J - t v + v 2 ) A™ 



A™*H™(5) with *?(S)=\n 



In 



1 + 



m 



5^B'Z- 2 / ? > 



Sfj,/m 2 + y/e 

£'=1429 for Z > 2 and £'=446 for Z=l. 

The maximum energy transferred to the photon is 

m(E — p) 
E(E — p + m) 

On the output all of the energy lost by a muon is assigned to the bremsstrahlung 
photon. 

9.3 Photonuclear interaction 




9.3.1 Bezrukov Bugaev parameterization of the photonuclear interaction 

The soft part of the photonuclear cross section is used as parametrized in [45] (un- 
derlined terms taken from [49,30] are important for tau propagation): 
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Fig. 37. Photon-nucleon cross sections, 
as described in the text: Kokoulin [43], 
W. Rhode [44], BB 1981 [45], ZEUS 94 
[46], ALLM 91 and 97 [47], Butkevich 
[48]. Curves 5-7 are calculated according 
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Fig. 38. Photonuclear energy losses (di- 
vided by energy), according to formulae 
from Section 9.3. Higher lines for the 
parameterizations 1-4 include the hard 
component [49], higher lines for 5-7 cal- 
culate shadowing effects as in Section 
9.3.3, lower as in Section 9.3.2 
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Nucleon shadowing is calculated according to 
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Several parametrization schemes for the photon-nucleon cross section are imple- 
mented. The first is 
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a- , (//) = 9G. L + -j=, for v < 17 GeV 



a lN {v) = 114.3 + 1.647 ln 2 [0.0213z/] fib, for v e [17, 200 GeV] [45] 
a lN (v) = 49.2 + 11.1 ln[z/] + 151. 8/ ^b, above 200 GeV [43]. 

The second is based on the table parametrization of [44] below 17 GeV. Since 
the second formula from above is valid for energies up to 10 6 GeV, it is taken to 
describe the whole energy range alone as the third case. Formula [46] 

a lN (v) = 63.5s - 097 + 145s" - 5 pb with s = 2Mv 

can also be used in the whole energy range, representing the fourth case (see Figure 
37). Finally, the ALLM parametrization (discussed in Section 9.3.2) or Butkevich- 
Mikhailov parameterization (discussed in Section 9.3.3) can be enabled. It does not 
rely on "nearly-real" exchange photon assumption and involves integration over 
the square of the photon 4-momentum (Q 2 ). Also, treatment of the hard component 
within the Bezrukov-Bugaev parameterization can optionally be enabled. The hard 
component of photonuclear cross section was calculated in [49] and parametrized 
in [30] as 



da h 



ard 



I 7 

= A--Y.dk logic v , used for 10_7 ^ v < 1. 1q2 GeV < E < 1q9 GeV - 



dv v k=Q 



Integration limits used for the photonuclear cross section are (kinematic limits for 
Q 2 are used for the ALLM and Butkevich-Mikhailov cross section formulae) 



mi M mi\ 

-^ w -^<Q 2 <2M(u-m^)-m 2 , E' = E-u. 



9.3.2 Abramowicz Levin Levy Maor (ALLM) parametrization of the photonuclear 
cross section 

The ALLM formula is based on the parametrization [50,47,51] 
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The limits of integration over Q 2 are given in the section for photonuclear cross 
section. 

F 2 = a(Z + (A-Z)P)F$ Here, a(A,x,Q 2 ) ~ a(A,x) 



a(A,x)=A-° A for x < 0.0014 

a(A, x) = A om9log ™ x+omr for 0.0014 < x < 0.04 

a{A,x) = l for x > 0.04 

P(x) = 1 - 1.85a; + 2.45a; 2 - 2.35a; 3 + x 4 
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(F? + F 2 R ) 



F*(x, Q 2 ) = ax? (1 - x) b > for i = P,R 



For / = c R , a R , b R , b P f(t) = f l + f 2 t h 

For g = c P , a P g(t) = g x + (g x - g 2 ) 



1 + 



i n Q 2 +Ql 
In A2 



In ^ 

111 A 2 



Q 2 + m 2 



for 



P,R, 



Q 2 + m 2 + W 2 -M 2 

where W is the invariant mass of the nucleus plus virtual photon [52]: W 2 = M 2 + 
2MEv — Q 2 . Figure 38 compares ALLM-parametrized cross section with formulae 
of Bezrukov and Bugaev from Section 9.3.1. 

The quantity R(x, Q 2 ) is not very well known, although it has been measured for 
high x (x > 0.1) [53] and modeled for small x (10~ 7 < x < 0.1, 0.01 GeV 2 < 
Q 2 < 50 GeV 2 ) [54]. It is of the order ~ 0.1 - 0.3 and even smaller for small Q 2 
(behaves as 0(Q 2 )). In Figure 39 three photonuclear energy loss curves for R=0, 
0.3, and 0.5 are shown. The difference between the curves never exceeds 7%. In 
the absence of a convenient parametrization for R at the moment, it is set to zero in 
MMC. 



The values of cross sections in Figures 37—39 should not be trusted at energies 
below 10 GeV. However, their exact values at these energies are not important for 
the muon propagation since the contribution of the photonuclear cross section to 
the muon energy losses in this energy range is negligible. 
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9.3.3 Butkevich-Mikhailov parametrization of the photonuclear cross section 

Following the parameterization of the proton (p) and neutron (n) structure functions 
according to the CKMT model [55,48], 

Fr(x, Q 2 ) = F p s > n (x, Q 2 ) + Ff$(x, Q 2 ) 



Ff (*, Q 2 ) = A s x~^\l - (^-) 
F2(x, Q 2 ) = A s x-^ 2 \1 - xT^ 

F P NS { X) Q 2 ) = xU v (x, Q 2 ) + xD v {x, Q 2 ) 
F£ s (x,Q 2 ) = - A xU v {x, Q 2 ) + AxD v (x, Q 2 ) 



1+A(Q 2 ) 



2 \ 1+A(Q 2 



xU v (x,Q 2 ) = B u x^- aR \\ - x) n{Q2) ^ 
xD v (x, Q 2 ) = B d x [1 - aR \l - x) n ^ +1 ^ 



2 \ a R 



Q 2 + bJ 

i ( Q 2 Y R 
Q 2 + b 



where A(Q 2 ) = A (l + -^-Y and n{Q 2 ) = \ ( 1 + Q ' 



Q 2 + d r y ^ J 2 \ Q 2 + c 



F A (x, Q 2 ) = r A ^ d [ZFi(x, Q 2 ) + (A - Z)F?(x, Q 2 )] 

The nuclear structure function r A l d can be evaluated as the shadowing function a 
from the previous section, or can optionally be calculated as follows [48]. At x > 

0.3, r A / d = 1 - m b (A)a osc (x), with m b (A) = M b [l - N S (A)/A\ and M b = 0.437. 
N S (A) is the Wood-Saxon potential 



N s {A)=Aix PQ J - 



r 2 dr 

1 + exp|f 

MA) 



+ exp{[r -r Q {A)]/aY 



where p = 0.17 far 3 , a = 0.54 fm, and r (A) = 1.12A 1 / 3 - 0.86A- 1 / 3 . 
a osc (x) = (1 - Xx) { (- - -) - /I ( — - \) ) , 



U CJ \u 2 c 2 



where u = 1 — x, c = 1 — x 2 , x 2 = 0.278, A = 0.5, and ji = m^/M. 
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At 10~ 3 > x < x ^ 0.3, r A / d (x) = x mi (l + m 2 )(l - m 3 x) with m, = Mjl - 
N S {A)/A\, where Mi = 0.129, M 2 = 0.456, and M 3 = 0.553. Here 



whereG(z/) is given by expression of Section 9.3.1 wither^ = 112.2(0.609z/ 00988 + 
1.037z/-°- 5944 ). At x < x function r A / d (x) = r A / d (x ). 



9.4 Electron pair production 



Two out of four diagrams describing pair production are shown below. These de- 
scribe the dominant "electron" term. The two diagrams not shown here describe the 
muon interacting with the atom and represent the "muon" term. The cross section 
formulae used here were first derived in [56,57,58]. 



1 



-i l/mi 



.1 + m 2 



(0.75G(u) + 0.25) 




Pt 




v = 



(e + + e_)/£, p=(e + -e_)/E 
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$ e = | [(2 + p 2 )(l + (3) + £(3 + p 2 )] In (l + f\ + 1 ^ ~ " (3 + p 2 )} ^ 



:i+p 2 ) (i + ^)-i(i+2/?)(i-p 2 ) 

+ g(1 ;^"^ ) +(i+2^)(i-p 2 )U M 



ln(l + 0+ 



^(i-p 2 ) 

3 m 



1 + 



2m v ^BZ- 1 / 3 (l+Q(l+y M ) 
£^(l-p 2 ) 

5 - p 2 + 4/3(1 + p 2 



2(1 + 3/3) ln(3 + 1/0 -P 2 - 2/3(2 - p 2 ) 

4 + p 2 + 3/3(1 + p 2 ) 

(l + p 2 )(3/2 + 2/3)ln(3 + + l-|p 2 



3m 
,2pT 



ZV3 ( i + e) (i + y e ; 



o _ ^ 2 t = (t lv \ 2 1 -P i 
P 2(1 -vY ^ \2m) 1-v 



0.073 lnf— - 

V 1+71 



Z 2 / 3 E/ fJ , 



0.26 



0.058 In ( 1+7 jfe g/ J - 0.14 



7i = 1.95 10~ 5 and 72 = 5.3 10~ 5 for Z^l 
7l = 4.4 10" 5 and 72 = 4.8 10~ 5 for Z = 1. 

Integration limits for this cross section are 



4m _ 3y/e p 

^ ^min — V — "Umax ^ A E* 



0< |p| < pn 



4m 



1 - 



Z l/3 

6p 2 



£ 2 (l-v) 



Muon pair production is discussed in detail in [59] and is not considered by MMC. 
Its cross section is estimated to be ~ 2 • 10 4 times smaller than the direct electron 
pair production cross section discussed above. 



9.5 Landau-Pomeranchuk-Migdal and Ter-Mikaelian effects 



These affect bremsstrahlung and pair production. See Figure 40 for the combined 
effect in ice and Frejus rock. 
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9.5.1 LPM suppression of the bremsstrahlung cross section: 

The bremsstrahlung cross section is modified as follows [60,61,62,63]: 

i(l_ w ) +T ,2_>M ^G(s) + 2[1 + (1 - vYMs)) . 

The regions of the following expressions for <f>(s) and G(s) were chosen to repre- 
sent the best continuous approximation to the actual functions: 



s 3 



^ 1 - 6XP ^ + ^ - ^ + 0.623 + 0.796, + 0.658,' ** s < 1-54954 



0(s) = l -0.012/s 4 for s> 1.54954 

/ 8s 2 
■?/>(s) = 1 — exp —4s 



1 + 3.936s + 4.97s 2 - 0.05s 3 + 7.50s 4 / 
G(s) = 3^(s) - 20(s) for s < 0.710390 
G(s) = 36s 2 /(36s 2 + 1) for 0.710390 < s < 0.904912 
G(s) = 1 - 0.022/s 4 for s > 0.904912. 

Here the SEB scheme [64] is employed for evaluation of 0(s), ip(s), and £(s) be- 
low: 



f(s') = 2 for s' <si 

j- / f\ 1 , 0.08(1 
£(s') = l + /i ^ 

£(s') = 1 for s' > 1 



, 0.08(1 - h)\l - (1 - h) 2 } r 

Ms') = l + h i ^ i for sx < s' < 1 

msi 



_ a(/zc 2 ) 2 X 

X is the same as in Section 9.7. Here are the rest of the definitions: 

s' r-Z l l z D n m e , I E LPM v Ins' 

s = —H: Si = V2 s — \ —— h = . 

B fj, ^8E{l-v) lnsi 



9.5.2 Dielectric (Longitudinal) suppression effect: 

In addition to the above change of the bremsstrahlung cross section, s is replaced 
by T • s and functions £(s), 0(s), and G(s) are scaled as [62] 

£(s)^£(rs) 0( s )^ (rs)/r g(s) -> G(rs)/r 2 . 
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1 ALLM R=0 

2 ALLM R=0.3 

3 ALLM R=0.5 
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ep.nr 




brems 



t'pair 



brems 
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Fig. 39. Comparison of ALLM en- 
ergy loss (divided by energy) for R=0 
(dashed-dotted), R=0.3 (dotted), #=0.5 
(dashed) 



Fig. 40. LPM effect in ice (higher plots) 
and Frejus rock (lower plots, multiplied 
by 10~ 3 ) 



Therefore the first formula in the previous section is modified as 

i( 1 -,) + ^^(^ + 2 [ i + (1 -^). 

T is defined as r = 1 + 7 2 (^^j , 

where uo p = JAitNZe 2 /m is the plasma frequency of the medium and vE is the 
photon energy. The dielectric suppression affects only processes with small photon 
transfer energy, therefore it is not directly applicable to the direct pair production 
suppression. 

9.5.3 LPM suppression of the direct pair production cross section: 
$ e from the pair production cross section is modified as follows [62,65]: 

$ e - ((1 +p)(A+[l + p 2 }B) + /3{C + [1 + p 2 ]D) + (1 - p 2 )E) ■ L e 



_ 1 \E L p M I 



4V v(l-p 2 Y 
The E LPM energy definition is different than in the brems strahlung case: 

E LPM = - ri fl 72r , where L = ln(3.25SZ -1 / 3 ). 
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Functions A(s,£), B(s,£), C(s, £), and D(s,£) are based on the approximation 
formulae 

\ 6s . m \ ( 6s ) 2 

$(s) = and G(s) 



6s + 1 w (6s) 2 + 1 

and are given below: 



G m r% s~i \ l 36s 2 (l + x) 2 + l ... 
A(s,x) = -(l + 2Gx)ln^^l G 



„_ / 36s 2 -1 

+6Gs 1 + ^T 



:j-x^ ^arctan(6s[x + 1]) — — 



5(s, x) = $(1 + $x) In 6s(1 + :r) + 1 - $ 

6sx 



„, , „ 2 . 36s 2 (l + x) 2 + l _ G 2 (36s 2 -1) / . . tT 

C(s, x) = -G 2 x In hG -x arctan(6s x + 1 ) 

36s z x z 6s V 



D(s,x) = $ - $ xln 



, (i.-i 1 • .rj • 1 

6sx 

7T 



E(s,x) = —6s ^arctan(6s[x + 1]) 



9.(5 Muon and tau decay 



Muon decay probability is calculated according to 

dN 1 



fix 7/3cr 

The energy of the outgoing electron is evaluated as 



V e = 7 ("rest + Purest ~ ™l cos (°) 



The value of cos(#) is distributed uniformly on (—1,1) and u rest is determined at 
random from the distribution 

^ = ^1(3-2^, X= JL. with ,„ in = m< md 
dx 1927T* z/ ma:r 2/x 



Tau leptonic decays, into a muon (17.37%) and electron (17.83%), are treated sim- 
ilarily. Hardronic decays are approximated by two-body decays into a neutrino 
and a hardonic part, which is assumed to be one of the particles or resonances: n 
(11.09%), p-770 (25.40%, M = 769.3 MeV), ai-1260 (18.26%, M = 1230 MeV), 
and the rest into p-1465 (10.05%, M = 1465 MeV). The energy of the hardronic 
part in the tau rest frame is evaluated as u rest = (m 2 + M 2 )/ (2m T ). 
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Fig. 41. Moliere scattering of one hundred 10 TeV muons going straight down 
through ice 

9. 7 Moliere scattering 



After passing through a distance x, the angular distribution is assumed Gaussian 
with a width y/26 [31,66,67]: 



#o = U - 6MeV zJx~/X Q [1 + 0.038 Hx/X Q )\ 
pep v 



X is evaluated as X 



n -l 



big) 



E 



big 



for E Ug ~ 10 20 eV 



Deviations in two directions perpendicular to the muon track are independent, but 
for each direction the exit angle and lateral deviation are correlated: 



Vpiane = z^Oq/v 12 + z 2 x9 /2 and 6 plane = z 2 9 




for independent standard Gaussian random variables (zi, z 2 ). A more precise treat- 
ment should take the finite size of the nucleus into account as described in [68]. See 
Figure 41 for an example of Moliere scattering of a high energy muon. 
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10 Conclusions 



A very versatile, clearly coded, and easy-to-use muon propagation Monte Carlo 
program (MMC) is presented. It is capable of propagating muon and tau leptons of 
energies from 105.7 MeV (muon rest mass, higher for tau) to 10 11 GeV (or higher), 
which should be sufficient for the use as propagator in the simulations of the mod- 
ern neutrino detectors. A very straightforward error control model is implemented, 
which results in computational errors being much smaller than uncertainties in the 
formulae used for evaluation of cross sections. It is very easy to "plug in" cross 
sections, modify them, or test their performance. The program was extended on 
many occasions to include new formulae or effects. MMC propagates particles in 
three dimensions and takes into account Moliere scattering on the atomic centers, 
which could be considered as the zeroth order approximation to true muon scat- 
tering since bremsstrahlung and pair production are effects that appear on top of 
such scattering. A more advanced angular dependence of the cross sections can be 
implemented at a later date, if necessary. 



Having been written in Java, MMC comes with the C/C++ interface package, which 
simplifies its integration into the simulation programs written in native languages. 
The distribution of MMC also includes a demonstration applet, which allows one 
to immediately visualize simulated events. 



MMC was incorporated into the simulation of the AMANDA, IceCube, and Frejus 
experiments. It is distributed at [9] with hope that the combination of precision, 
code clarity, speed, and stability will make it a useful tool in research where prop- 
agation of high energy particles through matter needs to be simulated. 



A calculation of coefficients in the energy loss formula dE/dx = a + bE and a 
similar formula for average range is presented for continuous (for energy loss) and 
stochastic (for average range calculation) energy loss treatments in ice and Frejus 
Rock. The calculated coefficients apply in the energy range from 20 GeV to 10 11 
GeV with an average deviation from the linear formula of 3.7% and maximum of 
6.6%. Also, 99.9% range of muons propagating in ice is estimated for energies 
from 1 GeV to 10 11 GeV. 



This work was supported by the German Academic Exchange Service (DAAD), 
U.S. NSF Grants OPP-020311 and OPP-0236449, and the U.S. DOE contract DE- 
AC-76SF00098. 
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A Tables used by Muon Monte Carlo (MMC) 



All cross sections were translated to units [1/cm] via multiplication by the number 
of molecules per unit volume. Many unit conversions (like eV — > J) were achieved 
using values of a = e 2 /he and r e = e 2 /m e c 2 . 



Summary of physical constants employed by MMC (as in [69]) 



a 1/137.03599976 


r e 2.817940285 • 10~ 13 cm 


N a 6.02214199 • 10 23 1/mol 


K 0.307075 MeV ■ cm 2 / g 


c 299792458 • 10 10 cm/s 


R y 13.60569172 eV 


m e 0.510998902 MeV 


m w 139.57018 Mev 


m p 938.271998 MeV 


m n 939.56533 MeV 


m M 105.658389 MeV 


2.19703- 10~ 6 s 


m T 1777.03 MeV 


r T 290.6 • 10~ 15 s 



Radiation logarithm constant B (taken from [70]) 



z 


B 


Z 


B 


Z 


B 


Z 


B 


Z 


B 


1 


202.4 


8 


\73.4 


15 


172.2 


22 


176.8 


53 


178.6 


2 


151.9 


9 


170.0 


16 


173.4 


26 


175.8 


74 


177.6 


3 


159.9 


10 


165.8 


17 


174.3 


29 


173.1 


82 


178.0 


4 


172.3 


11 


165.8 


18 


174.8 


32 


173.0 


92 


179.8 


5 


177.9 


12 


167.1 


19 


175.1 


35 


173.5 






6 


178.3 


13 


169.1 


20 


175.6 


42 


175.9 


other 


182.7 


7 


176.6 


14 


170.8 


21 


176.2 


50 


177.4 
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Media constants (taken from [6,33]) 



Material 


Z 


A 


I, eV 


-C 


a 


m 


^0 


Xi 


p, g/cm 2 


So 


Water 


1 + 


1.00794 


75.0 


3.5017 


0.09116 


3.4773 


0.2400 


2.8004 


1.000 





Ice 


+ 8 


15.9994 


75.0 


3.5017 


0.09116 


3.4773 


0.2400 


2.8004 


0.917 





Stand. Rock 


11 


22 


136.4 


3.7738 


0.08301 


3.4120 


0.0492 


3.0549 


2.650 





Frejus Rock 


10.12 


20.34 


149.0 


5.053 


0.078 


3.645 


0.288 


3.196 


2.740 





Iron 


26 


55.845 


286.0 


4.2911 


0.14680 


2.9632 


-0.0012 


3.1531 


7.874 


0.12 


Hydrogen 


1 


1.00794 


21.8 


3.0977 


0.13483 


5.6249 


0.4400 


1.8856 


0.07080 





Lead 


82 


207.2 


823.0 


6.2018 


0.09359 


3.1608 


0.3776 


3.8073 


11.350 


0.14 


Uranium 


92 


238.0289 


890.0 


5.8694 


0.19677 


2.8171 


0.2260 


3.3721 


18.950 


0.14 


N 2 78.1% 7 


14.0067 


















< 2 21.0% 8 


15.9994 


85.7 


10.5961 


0.10914 


3.3994 


1.7418 


4.2759 


1.205 

•10~ 3 





Ar 0.9% 18 


39.948 

















Parameterization coefficients of the hard component of the photonuclear cross section (as in [30]) 



E 


10 3 GeV 


10 4 GeV 


10 5 GeV 


10 6 GeV 


10 7 GeV 


10 8 GeV 


10 9 GeV 


muons 


a 


7.174409 • 10" 4 


1.7132 • 10~ 3 


4.082304 • 10~ 3 


8.628455 • 10~ 3 


0.01244159 


0.02204591 


0.03228755 


ai 


-0.2436045 


-0.5756682 


-1.553973 


-3.251305 


-5.976818 


-9.495636 


-13.92918 


a 2 


-0.2942209 


-0.68615 


-2.004218 


-3.999623 


-6.855045 


-10.05705 


-14.37232 


a 3 


-0.1658391 


-0.3825223 


-1.207777 


-2.33175 


-3.88775 


-5.636636 


-8.418409 


04 


-0.05227727 


-0.1196482 


-0.4033373 


-0.7614046 


-1.270677 


-1.883845 


-2.948277 


a 5 


-9.328318 • IO- 3 


-0.02124577 


-0.07555636 


-0.1402496 


-0.2370768 


-0.3614146 


-0.5819409 


a 6 


-8.751909 • 10~ 4 


-1.987841 • 10~ 3 


-7.399682 • 10~ 3 


-0.01354059 


-0.02325118 


-0.03629659 


-0.059275 


a 7 


-3.343145 • 10~ 5 


-7.584046 • 10~ 5 


-2.943396 • 10~ 4 


-5.3155 • 10~ 4 


-9.265136 • 10~ 4 


-1.473118 • 10~ 3 


-2.419946 • 10~ 3 


taus 


a 


-1.269205 • 10~ 4 


-2.843877 • 10" 4 


-5.761546 • 10~ 4 


-1.195445 • 10~ 3 


-1.317386 • 10~ 3 


-9.689228 • 10~ 15 


-6.4595 • 10~ 15 


ai 


-0.01563032 


-0.03589573 


-0.07768545 


-0.157375 


-0.2720009 


-0.4186136 


-0.8045046 




0.04693954 


0.1162945 


0.3064255 


0.7041273 


1.440518 


2.533355 


3.217832 


a 3 


0.05338546 


0.130975 


0.3410341 


0.7529364 


1.425927 


2.284968 


2.5487 


04 


0.02240132 


0.05496 


0.144945 


0.3119032 


0.5576727 


0.8360727 


0.8085682 


a 5 


4.658909 • 10~ 3 


0.01146659 


0.03090286 


0.06514455 


0.1109868 


0.1589677 


0.1344223 


a 6 


4.822364 • 10~ 4 


1.193018 • 10~ 3 


3.302773 • 10~ 3 


6.843364 • 10~ 3 


0.011191 


0.015614 


0.01173827 


a 7 


1.9837 • 10~ 5 


4.940182 • 10~ 5 


1.409573 • 10~ 4 


2.877909 • 10~ 4 


4.544877 • 10~ 4 


6.280818 • 10~ 4 


4.281932 • 10~ 4 



ALLM ('91) parameters (as in [50,71]) 



a pi 


-0.04503 


&P2 


-0.36407 


ap3 


8.17091 


am 


0.60408 




0.17353 


dR3 


1.61812 


bpi 


0.49222 2 


bp 2 


0.52116 2 


bp 3 


3.55115 


bm 


1.26066 2 


bR2 


1.83624 2 


bR3 


0.81141 


cpi 


0.26550 


CP2 


0.04856 


CP3 


1.04682 


cm 


0.67639 


CR2 


0.49027 


CR3 


2.66275 


mp 


10.67564 • 10 6 MeV 2 


A 2 


0.06527 • 10 6 MeV 2 


m\ 


0.30508 • 10 6 MeV 2 


m\ 


0.20623 • 10 6 MeV 2 


Ql - A 2 


0.27799 • 10 6 MeV 2 







ALLM ('97) parameters (as in [47,71]) 



a pi 


-0.0808 


op 2 


-0.44812 


ap3 


1.1709 


ClRl 


0.58400 




0.37888 


dR3 


2.6063 


bpi 


0.60243 2 


bp2 


1.3754 2 


bp 3 


1.8439 


bm 


0.10711 2 


b R 2 


1.9386 2 


bp3 


0.49338 


Cpi 


0.28067 


CP2 


0.22291 


CP3 


2.1979 


CRl 


0.80107 


CR2 


0.97307 


CR3 


3.4942 


mp 


49.457 • 10 6 MeV 2 


A 2 


0.06527- 10 6 MeV 2 


ml 


0.31985 • 10 6 MeV 2 


m\ 


0.15052 • 10 6 MeV 2 


Ql - A 2 


0.46017- 10 6 MeV 2 







CKMT parameters of the Butkevich-Mikhailov parameterization (as in [48,72]) 



a 


0.2513 


• 10 6 


MeV 2 


A 


0.0988 


b 


0.6186 


• 10 6 


MeV 2 


CiR 


0.4056 


c 


3.0292 


• 10 6 


MeV 2 


T 


1.8152 


d 


1.4817 


• 10 6 


MeV 2 


B u 


1.2437 






0.12 




B d 


0.1853 
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Tsai's Radiation logarithms L ra d and L' rad (as in [31,40]) 



z 


L ra d 


n rad 


1 


5.31 


6.144 


2 


4.79 


5.621 


3 


4.74 


5.805 


4 


4.71 


5.924 


other 


ln(184.15Z- 1 / 3 ) 


ln(1194Z~ 2 / 3 ) 



B Comparison of Spectra of Secondaries Produced with MMC, 
MUM [12], LOH [6], and LIP [7,17] 



Fig. B.l. Spectra of the secondaries: the setup 

Az 



20 m 



B.l Spectra of the secondaries 



In order to determine spectra of primaries consistently for all programs, the fol- 
lowing setup was used. For each muon with fixed initial energy a first secondary 
created within the first 20 meters is recorded (Figure B.l). This is somewhat differ- 
ent from what was done for Figure 7, since the energy of the muon at the moment 
when the secondary is created is somewhat smaller than the initial energy due to 
continuous energy losses. These are smaller when v cut is smaller, and are generally 
negligible for all cases considered below. 

In Figure B.2 solid curves are probability functions normalized to the total number 
of secondaries above 500 MeV. In Figure B.4 solid curves are probability func- 
tions normalized to the total number of secondaries above 1CT 3 ■ E^. In Figure B.5 



41 



2 3 2 3 4 5 6 7 8 9 2 4 6 S 10 12 14 16 16 20; 

1 10 10 10 1 10 10 10 10 10 10 10 10 10 1 10 10 10 10 10 10 10 10 10 1010 

energy [GeV] energy [GeV] energy [GeV] 



Fig. B.2. MMC: E^ = 10 3 GeV, = 10 9 GeV, and = 10 21 GeV with 
E cut = 500 MeV 
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Fig. B.3. MUM: = 10 3 GeV, = 10 6 GeV, and E^ = 10 9 GeV with 

E cut = 500 MeV 




Fig. B.4. LOH: E^ = 10 3 GeV, E^ = IO 4 GeV, and = 10 6 GeV 

solid curves are probability functions normalized to the total number of secondaries 
above 1CT 2 • E^. A setting of E big = 10 21 GeV is used for the third plot in Figure 
B.2 (default is 10 11 GeV). 



B.2 Number and total energy of secondaries 



In spite of the numerous problems with propagation codes other than MMC, shown 
in Figures B.3— B.5, it was possible to use these codes in the simulation of AMANDA- 
II. To understand why, the following setup is used. For each muon with fixed initial 
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Fig. B. 6. Number and total energy of sec- „ _ _ . T , . , . 

, . Fig. B.7. Number 01 secondaries 

ondanes: the setup 

energy all secondaries created within the first 800 meters (equal to the height of 
the AMANDA-II detector) are recorded (Figure B.6). Although the number of sec- 
ondaries generated by propagators LOH and LIP is different from that generated 
by MMC or MUM (Figure B.7), the total energy deposited in the volume of the 
detector is commensurable between all four propagators. The number of generated 
secondaries depends on the chosen value of E cut or v cut . While MMC and MUM 
allow one to select this value, LOH and LIP have a built-in value which cannot be 
changed. From Figure B.7 it appears that these codes use a value of v cut which lies 
between 1CT 2 and 10~ 3 since their number of secondaries lies between that gen- 
erated with MMC with v cut = 1CT 2 and v cut = 1CT 3 . One would expect the total 
energy of secondaries generated with LOH or LIP to be somewhat lower than that 
generated with MMC or MUM with E cut = 500 MeV. This, however, is not true: 
the total energy of secondaries generated with LOH and LIP is somehow renormal- 
ized to match that of MMC and MUM (Figures B.8 and B.9). 

Figures B.7 and B.9 also demonstrate the span of energies over which MMC can be 
used with fixed E cut = 0.5 GeV. With such E cut , MMC seems to work for energies 
up to 0.5 • 10 15 GeV, which is determined by the computer precision with which 
double precision numbers can be added: 0.5/0.5 • 10 15 ~ 10 -15 . When relative po- 
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Fig. B.8. Total energy of secondaries Fig. B.9. Relative energy of secondaries 

sition increments fall below that, the muon "gets stuck" in one point until its energy 
becomes sufficiently low or it propagates without stochastic losses sufficiently far, 
so that it can advance again. A muon "stuck" in this fashion still looses the energy, 
which is why it appears that its losses go up. With fixed v cut = 1CT 2 — 10~ 3 (and 
apparently as low as 1CT 12 — 1CT 15 ), MMC shows no signs of such deterioration. 
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